import os
import numpy as np
import pandas as pd
import nibabel as nib
from scipy.spatial.distance import cdist
from sklearn.metrics import pairwise_distances as pdist
import networkx as nx
import plotly.graph_objs as go
import plotly.figure_factory as ff
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import linecache
from itertools import compress
import xml.etree.ElementTree as ET
import matplotlib as mpl
from matplotlib import cm
# Color maps
viridis = cm.get_cmap('viridis')
tab10 = cm.get_cmap('tab10')
## To save as vector image use Orca, e.g.:
# /home/rhaast/.local/bin/orca graph vascular_tree_course.json -o vascular_tree_course -f pdf
Some global plotly settings
# Adjust camera position
camera_vasculature = dict(
eye=dict(x=0, y=0, z=1.5)
)
camera_hippocampus = dict(
eye=dict(x=0, y=0, z=2.5)
)
# Hide axis planes
scene = dict(
xaxis = dict(showbackground=False, visible=False),
yaxis = dict(showbackground=False, visible=False),
zaxis = dict(showbackground=False, visible=False)
)
# Set margins
margin = dict(l=0, r=0, b=0, t=0)
# Combine into master dictionary
vasculature_settings = dict(
width=675,
height=500,
scene=scene,
scene_aspectmode='data',
scene_camera=camera_vasculature,
margin=margin,
showlegend=False,
hovermode=False
)
hippocampus_settings = dict(
width=675,
height=500,
scene=scene,
scene_aspectmode='data',
scene_camera=camera_hippocampus,
margin=margin,
showlegend=False,
hovermode=False
)
subject = os.environ["SUBJECT"]
hemisphere = os.environ["HEMISPHERE"]
# subject = '01'
# hemisphere = 'R'
mri_maps = ['B1map','PWI']
vasculature_path = '../results/vasculature/sub-{}/{}'.format(subject, hemisphere)
in_xml = '{}/vascular_tree.xml'.format(vasculature_path)
in_centerline = '{}/vascular_tree_centerline.nii.gz'.format(vasculature_path)
in_xfm = '../resources/mevislabSTL2Orig.txt'
in_lores_surface = '../results/surface_warps/sub-{}/{}/midthickness.downsampled.native.surf.gii'.format(subject, hemisphere)
in_hires_surface = '../results/surface_warps/sub-{}/{}/midthickness.native.surf.gii'.format(subject, hemisphere)
mri_path = '../results/maps/sub-{}'.format(subject, hemisphere)
in_mri = '{}/sub-{}_{}_{}.nii.gz'
parameters = linecache.getline(in_xfm, 4)
parameters = parameters.split(' ')[1:10]
xfm = np.array([ float(p) for p in parameters ]).reshape(3,3)
tree = ET.parse(in_xml)
root = tree.getroot()
# Lores surface, for visualization
surf = nib.load(in_lores_surface)
vertices_lores = surf.darrays[0].data
faces_lores = surf.darrays[1].data
n_lores_vertices = len(vertices_lores)
# Hires surface, for output maps
surf = nib.load(in_hires_surface)
vertices_hires = surf.darrays[1].data
faces_hires = surf.darrays[0].data
n_hires_vertices = len(vertices_hires)
nii = nib.load(in_centerline)
centerline = nii.get_fdata()
xfm_nii = nii.affine
voxel_size = abs(nii.affine[0,0])
centerline_idx = np.argwhere(centerline)
centerline_coords = [ np.dot(np.hstack((centerline_idx[i],1)),xfm_nii.T)[:3] for i in range(0,len(centerline_idx)) ]
centerline_coords = np.array(centerline_coords)
mri_data = np.zeros((centerline.shape[0],centerline.shape[1],centerline.shape[2],len(mri_maps)))
for p,parameter in enumerate(mri_maps):
in_nii = in_mri.format(mri_path, subject, parameter, hemisphere)
nii_data = nib.load(in_nii).get_fdata()
mri_data[1:-1,1:-1,1:-1,p] = nii_data
x = []
y = []
z = []
d = []
for skeleton in root.iter('Skeleton'):
coords = skeleton.find('_PropertyContainer').find('PropertyValues').find('Vector3').text
x.append(float(coords.split(' ')[0]))
y.append(float(coords.split(' ')[1]))
z.append(float(coords.split(' ')[2]))
diameter = skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text
d.append(float(diameter)*2)
data_skeleton = np.vstack((x,y,z,d))
data_skeleton[:3,:] = np.dot(xfm,data_skeleton[:3,:])
edge_diameters = []
for skeleton in root.iter('Skeletons'):
d = []
for _skeleton in skeleton.iter('Skeleton'):
diameter = _skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text
d.append(float(diameter))
edge_diameters.append(np.median(d))
node_id = 0
df_main_nodes = pd.DataFrame(columns=['x','y','z','diameter'])
for node in root.iter('node'):
coords = node.find('_BaseGraphItem').find('_PropertyContainer').find('PropertyValues').find('Vector3').text
node_x = float(coords.split(' ')[0])
node_y = float(coords.split(' ')[1])
node_z = float(coords.split(' ')[2])
node_xfm = np.dot(xfm,np.array((node_x,node_y,node_z)))
diameter = node.find('_BaseGraphItem').find('_PropertyContainer').find('PropertyValues').find('double').text
df_main_nodes.loc[node_id] = [node_xfm[0], node_xfm[1], node_xfm[2], float(diameter)*2]
node_id += 1
Only contains edges between the main (i.e., starting, root and branching) nodes.
# MeVisLab uses 1-based index, so minus one
pre_node = [int(edge.find('PredId').text)-1 for edge in root.iter('edge')]
post_node = [int(edge.find('SuccId').text)-1 for edge in root.iter('edge')]
df_main_edges = pd.DataFrame(list(zip(pre_node,post_node)), columns=['pre','post'])
Uses skeleton elements from XML data
# Initate pandas dataframe for nodes and edges
node_id = 0
df_nodes = pd.DataFrame(columns=['x','y','z','diameter','b1','pwi'])
edge_id = 0
df_edges = pd.DataFrame(columns=['pre','post','distance'])
# Itererate through edges
for edge in root.iter('edge'):
# Iterate through elements of skeleton
ids = None
per_skeleton = None
for skeleton in edge.iter('Skeleton'):
# Parse x,y,z coordinates
skel_xyz = skeleton.find('_PropertyContainer').find('PropertyValues').find('Vector3').text
skel_x = float(skel_xyz.split(' ')[0])
skel_y = float(skel_xyz.split(' ')[1])
skel_z = float(skel_xyz.split(' ')[2])
# Transform coordinates
skel_xfm = np.dot(xfm,np.array((skel_x,skel_y,skel_z)))
# Parse diameter
diameter = skeleton.find('_PropertyContainer').find('PropertyValues').find('double').text
# Extract MRI-based data
parameter_data = []
for p,parameter in enumerate(mri_maps):
tmp = np.sum(np.abs(centerline_coords - [skel_xfm[0], skel_xfm[1], skel_xfm[2]]),axis=1)
idx = centerline_idx[np.argmin(tmp)]
parameter_data.append(mri_data[idx[0],idx[1],idx[2],p])
df_nodes.loc[node_id] = [skel_xfm[0], skel_xfm[1], skel_xfm[2], float(diameter)*2, parameter_data[0], parameter_data[1]]
# Increase id
ids = np.hstack((ids, node_id)) if ids is not None else node_id
node_id += 1
# Define edges based on collected id's
for i in ids[:-1]:
distance = cdist(
np.array((df_nodes.x[i],df_nodes.y[i],df_nodes.z[i])).reshape(-1,1),
np.array((df_nodes.x[i+1],df_nodes.y[i+1],df_nodes.z[i+1])).reshape(-1,1)
)
df_edges.loc[edge_id] = [i,i+1,distance[0][0]]
edge_id += 1
Dataframes need to be cleaned as there will be duplicates near the main nodes
duplicate_entries = df_nodes[df_nodes.duplicated(['x','y','z','diameter'])]
for i in range(0,len(duplicate_entries)):
duplicates = df_nodes[(df_nodes.x == duplicate_entries.iloc[i].x) &
(df_nodes.y == duplicate_entries.iloc[i].y) &
(df_nodes.z == duplicate_entries.iloc[i].z) &
(df_nodes.diameter == duplicate_entries.iloc[i].diameter)
]
for d in duplicates.index.values[1:]:
df_edges = df_edges.replace(d, duplicates.index.values[0])
df_nodes = df_nodes.drop(d)
G_nodes = df_main_nodes.index.values
G_course = nx.Graph()
G_course.add_nodes_from(G_nodes)
for n in G_nodes:
G_course.nodes[n]['position'] = (df_main_nodes.x[n], df_main_nodes.y[n], df_main_nodes.z[n])
G_course.nodes[n]['diameter'] = (df_main_nodes.diameter[n])
for e in np.arange(0,df_main_edges.shape[0]):
G_course.add_edge(
int(df_main_edges.iloc[e]['pre']),
int(df_main_edges.iloc[e]['post'])
)
nodes = np.unique(df_main_nodes.index.values)
G_course_props_dict = {n: (
df_main_nodes.x[n], df_main_nodes.y[n], df_main_nodes.z[n], df_main_nodes.diameter[n], G_course.degree[n]
) for n in nodes}
node_degree_color_dict = {1 : 'green', 2 : 'black', 3 : 'red', 4 : 'purple', 5 : 'blue'}
node_degree_color = [None] * nx.number_of_nodes(G_course)
for node, degree in G_course.degree:
node_degree_color[node] = node_degree_color_dict[degree]
fig = go.Figure()
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
color='lightgray')])
fig.add_trace(go.Scatter3d(
x=df_main_nodes.x,
y=df_main_nodes.y,
z=df_main_nodes.z,
mode='markers',
marker=dict(
size=3,
color=node_degree_color, #>'black'
)))
for i,j in enumerate(G_course.edges()):
x = np.array((G_course_props_dict[j[0]][0], G_course_props_dict[j[1]][0]))
y = np.array((G_course_props_dict[j[0]][1], G_course_props_dict[j[1]][1]))
z = np.array((G_course_props_dict[j[0]][2], G_course_props_dict[j[1]][2]))
fig.add_trace(
go.Scatter3d(
x=x,
y=y,
z=z,
mode='lines',
line=dict(
width=3,
color='black'
)
)
)
fig.update_layout(dict1=vasculature_settings)
fig.write_json("{}/vascular_tree_course.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_course.html".format(vasculature_path))
fig.show()
Used cleaned dataframes to define network
G_nodes = df_nodes.index.values
G_fine = nx.Graph()
G_fine.add_nodes_from(G_nodes)
for n in G_nodes:
G_fine.nodes[n]['position'] = (df_nodes.x[n], df_nodes.y[n], df_nodes.z[n])
G_fine.nodes[n]['diameter'] = (df_nodes.diameter[n])
for e in np.arange(0,df_edges.shape[0]):
G_fine.add_edge(
int(df_edges.iloc[e]['pre']),
int(df_edges.iloc[e]['post']),
distance=float(df_edges.iloc[e]['distance'])
)
nodes = np.unique(df_nodes.index.values)
G_fine_props_dict = {n: (
df_nodes.x[n], df_nodes.y[n], df_nodes.z[n], df_nodes.diameter[n], G_fine.degree[n]
) for n in nodes}
fig = go.Figure()
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
color='lightgray')])
fig.add_trace(go.Scatter3d(
x=df_main_nodes.x,
y=df_main_nodes.y,
z=df_main_nodes.z,
mode='markers',
marker=dict(
size=3,
color=node_degree_color, #>'black'
)))
for i,j in enumerate(G_fine.edges()):
x = np.array((G_fine_props_dict[j[0]][0], G_fine_props_dict[j[1]][0]))
y = np.array((G_fine_props_dict[j[0]][1], G_fine_props_dict[j[1]][1]))
z = np.array((G_fine_props_dict[j[0]][2], G_fine_props_dict[j[1]][2]))
fig.add_trace(
go.Scatter3d(
x=x,
y=y,
z=z,
mode='lines',
line=dict(
width=10,
color='black'
)
)
)
fig.update_layout(dict1=vasculature_settings)
fig.write_json("{}/vascular_tree_fine.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_fine.html".format(vasculature_path))
fig.show()
Edges are treated as truncated cones with varying radii and height (i.e., radii and distance between adjacent nodes)
def calculate_edge_geometry(r1, r2, h):
import math
v = (1/3)*math.pi*(pow(r1,2)+r1*r2+pow(r2,2))*h # volume
la = math.pi*(r1+r2)*np.sqrt(pow((r1-r2),2)+pow(h,2)) # lateral area
sa = la + math.pi*(pow(r1,2)+(pow(r2,2))) # surface area
return(v, la, sa)
edge_geometry = [calculate_edge_geometry(
df_nodes.diameter[df_edges.pre[i]]/2, df_nodes.diameter[df_edges.pre[i]]/2, df_edges.distance[i]
) for i in df_edges.index.values]
edge_geometry = np.array(edge_geometry)
df_edges['volume'] = edge_geometry[:,0]
df_edges['lateral_area'] = edge_geometry[:,1]
df_edges['surface_area'] = edge_geometry[:,2]
# Project edge geometry data on neighboring nodes
node_geometry = np.zeros((nx.number_of_nodes(G_fine),6))
for n, node in enumerate(df_nodes.index.values):
df_filtered = df_edges[(df_edges.pre==node) | (df_edges.post==node)]
# Vessel geometry
node_geometry[n,0] = np.mean(df_filtered[df_filtered.volume!=0].volume)
node_geometry[n,1] = np.mean(df_filtered[df_filtered.volume!=0].lateral_area)
node_geometry[n,2] = np.mean(df_filtered.surface_area)
node_geometry[n,3] = np.mean(
df_nodes.loc[pd.unique(
df_filtered[['pre','post']].values.ravel()
)].diameter
)
# Smooth MRI-sampled data
node_geometry[n,4] = np.mean(
df_nodes.loc[pd.unique(
df_filtered[['pre','post']].values.ravel()
)].b1
)
node_geometry[n,5] = np.mean(
df_nodes.loc[pd.unique(
df_filtered[['pre','post']].values.ravel()
)].pwi
)
# Transform to hex codes for plotting
volume_norm = mpl.colors.Normalize(vmin=min(node_geometry[:,0]), vmax=max(node_geometry[:,0]))
volume_color = [ mpl.colors.to_hex(viridis(volume_norm(node_geometry[i,0]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]
lateral_area_norm = mpl.colors.Normalize(vmin=min(node_geometry[:,1]), vmax=max(node_geometry[:,1]))
lateral_area_color = [ mpl.colors.to_hex(viridis(lateral_area_norm(node_geometry[i,1]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]
surface_area_norm = mpl.colors.Normalize(vmin=min(node_geometry[:,2]), vmax=max(node_geometry[:,2]))
surface_area_color = [ mpl.colors.to_hex(viridis(surface_area_norm(node_geometry[i,2]))) for i in np.arange(0,nx.number_of_nodes(G_fine)) ]
from sklearn.preprocessing import minmax_scale
diameter_scaled = minmax_scale(node_geometry[:,3], feature_range=(3,10))
Identify connected components in graph and sort based on size
# Identify connected components and sort from largest to smallest (i.e., connected nodes)
G_cc = sorted(nx.connected_components(G_fine), key=len, reverse=True)
# Treshold for masking small components
treshold = 200
components_mask = np.ones((nx.number_of_nodes(G_fine)))
# Create compononts color map
components_color = [None] * nx.number_of_nodes(G_fine)
for c, component in enumerate(G_cc):
G_component = G_fine.subgraph(G_cc[c])
intersect, ind_a, ind_b = np.intersect1d(
list(G_fine.nodes), list(G_component.nodes), return_indices=True)
for i in ind_a:
components_color[i] = mpl.colors.to_hex(tab10(c))
# Mask if components smaller than treshold
if len(ind_a) <= 100:
components_mask[ind_a] = 0
components_mask = np.array(components_mask,dtype='bool')
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
color='lightgray',
showscale=False)])
fig.add_trace(go.Scatter3d(
x=df_nodes.x[components_mask],
y=df_nodes.y[components_mask],
z=df_nodes.z[components_mask],
mode='markers',
marker=dict(
size=diameter_scaled[components_mask],
line_width=0,
opacity=1,
color=node_geometry[components_mask,3], #list(compress(node_geometry[:,3], components_mask))
colorscale='plasma',
showscale=True,
colorbar=dict(
title="Diameter (mm)",
titleside="right",
thickness=10
)
)))
fig.update_layout(dict1=vasculature_settings)
fig.write_json("{}/vascular_tree_diameter.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_diameter.html".format(vasculature_path))
fig.show()
Generate shortest path from root nodes and export to matrix for color mapping
# Extract root nodes from XML
root_nodes = [int(r.text) for r in root.iter('RootID')]
# Identify root nodes based on x, y, z coordinates
root_node_ids = []
for root_node in root_nodes:
x = df_main_nodes.iloc[root_node-1].x
y = df_main_nodes.iloc[root_node-1].y
z = df_main_nodes.iloc[root_node-1].z
# x = df_main_nodes[1,df_main_nodes[0,:]==root_node][0]
# y = df_main_nodes[2,df_main_nodes[0,:]==root_node][0]
# z = df_main_nodes[3,df_main_nodes[0,:]==root_node][0]
df_root = df_nodes.query('x=={} & y=={} & z=={}'.format(x,y,z))
root_node_ids.append(df_root.index.values[0])
# Initiate array for filling with calculated distances
distance_to_root = np.zeros((len(G_fine.nodes),len(root_node_ids)))
for r, root_node in enumerate(root_node_ids):
# Define source node
source_node = root_node
# Calculate shortest path from source node to each other connected node
shortest_path_from_root = nx.shortest_path(G_fine, source=source_node, weight='distance')
# Iterate through each connected node and its path to the source node
for node, path in shortest_path_from_root.items():
# Find index of node
node_index = np.where(G_nodes == node)[0][0]
# Skip source node
if node != source_node:
mask = list(zip(path[:-1], path[1:]))
df_x = pd.MultiIndex.from_frame(df_edges[['pre','post']])
df_y = pd.MultiIndex.from_frame(df_edges[['post','pre']])
dist_x = sum(df_edges[df_x.isin(mask)].distance.values)
dist_y = sum(df_edges[df_y.isin(mask)].distance.values)
distance_to_root[node_index,r] = dist_x + dist_y
Most anterior, inferior vessel segment characterized by largest diameter
start_nodes = np.zeros((len(G_cc)))
for c, component in enumerate(G_cc):
order_x = df_nodes['x'][list(G_cc[c])].argsort()
rank_x = order_x.argsort()
order_y = df_nodes['y'][list(G_cc[c])].argsort()[::-1]
rank_y = order_y.argsort()
order_d = df_nodes['diameter'][list(G_cc[c])].argsort()[::-1]
rank_d = order_d.argsort()
# x, y, and diameter are weighted differentely.
rank_sum = np.average(np.vstack((rank_x,rank_y,rank_d)).T,weights=[0.5,1,1], axis=1)
start_nodes[c] = df_nodes.loc[list(G_cc[c])].iloc[np.argmin(rank_sum)].name
distance_to_start = np.zeros((len(G_fine.nodes),len(G_cc)))
for s, start_node in enumerate(start_nodes):
# Define source node
source_node = start_node
# Calculate shortest path from source node to each other connected node
shortest_path_from_root = nx.shortest_path(G_fine, source=source_node, weight='distance')
# Iterate through each connected node and its path to the source node
for node, path in shortest_path_from_root.items():
# Find index of node
node_index = np.where(G_nodes == node)[0][0]
# Skip source node
if node != source_node:
mask = list(zip(path[:-1], path[1:]))
df_x = pd.MultiIndex.from_frame(df_edges[['pre','post']])
df_y = pd.MultiIndex.from_frame(df_edges[['post','pre']])
dist_x = sum(df_edges[df_x.isin(mask)].distance.values)
dist_y = sum(df_edges[df_y.isin(mask)].distance.values)
distance_to_start[node_index,s] = dist_x + dist_y
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
color='lightgray',
showscale=False)])
fig.add_trace(go.Scatter3d(
x=df_nodes.x[components_mask],
y=df_nodes.y[components_mask],
z=df_nodes.z[components_mask],
mode='markers',
marker=dict(
size=diameter_scaled[components_mask],
line_width=0,
opacity=1,
color=np.max(distance_to_start[components_mask],axis=1),
colorscale='Viridis',
cmax=60,
cmin=0,
showscale=True,
colorbar=dict(
title="Geodesic distance from seed (mm, per component)",
titleside="right",
thickness=10
)
),
))
fig.update_layout(
dict1 = vasculature_settings,
scene=dict(
xaxis = dict(showbackground=False, visible=False),
yaxis = dict(showbackground=False, visible=False),
zaxis = dict(showbackground=False, visible=False),
annotations=[
dict(
x=df_nodes.loc[start_nodes[0]]['x'],
y=df_nodes.loc[start_nodes[0]]['y'],
z=df_nodes.loc[start_nodes[0]]['z'],
text="Seed",
bgcolor="black",
font=dict(
color="white",
size=16
),
showarrow=True,
arrowcolor="black"
)
])
)
fig.write_json("{}/vascular_tree_distance_to_seed.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_distance_to_seed.html".format(vasculature_path))
fig.show()
hippocampus_vessel_metrics_hires = np.zeros((n_hires_vertices,3))
hippocampus_vessel_metrics = np.zeros((n_lores_vertices,3))
# Hires for exporting
for i in range(0,n_hires_vertices):
dist = cdist(vertices_hires[i].reshape(1,3), df_nodes.values[:,:3])
hippocampus_vessel_metrics_hires[i,0] = df_nodes.values[:,3][np.argmin(dist)]
hippocampus_vessel_metrics_hires[i,1] = np.max(distance_to_start,axis=1)[np.argmin(dist)]
hippocampus_vessel_metrics_hires[i,2] = np.amin(dist)-(0.5*hippocampus_vessel_metrics_hires[i,0])
# Lores for exporting
for i in range(0,n_lores_vertices):
dist = cdist(vertices_lores[i].reshape(1,3), df_nodes.values[:,:3])
hippocampus_vessel_metrics[i,0] = df_nodes.values[:,3][np.argmin(dist)]
hippocampus_vessel_metrics[i,1] = np.max(distance_to_start,axis=1)[np.argmin(dist)]
hippocampus_vessel_metrics[i,2] = np.amin(dist)-(0.5*hippocampus_vessel_metrics[i,0])
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
intensity=hippocampus_vessel_metrics[:,2],
colorscale='magma',
showscale=True,
colorbar=dict(
title="Euclidean distance to closest vessel (mm)",
titleside="right",
thickness=10
)
)])
fig.update_layout(dict1=hippocampus_settings)
fig.write_json("{}/hippocampus_distance_to_vessel.json".format(vasculature_path))
fig.write_html("{}/hippocampus_distance_to_vessel.html".format(vasculature_path))
fig.show()
hippocampus_mri_metrics_hires = np.zeros((n_hires_vertices,len(mri_maps)))
hippocampus_mri_metrics = np.zeros((n_lores_vertices,len(mri_maps)))
# Hires for exporting
for i in range(0,n_hires_vertices):
idx = np.round(np.dot(np.hstack((vertices_hires[i,:],1)),np.linalg.inv(xfm_nii).T))
for p,parameter in enumerate(mri_maps):
hippocampus_mri_metrics_hires[i,p] = mri_data[int(idx[0]),int(idx[1]),int(idx[2]),p]
# Lores for exporting
for i in range(0,n_lores_vertices):
idx = np.round(np.dot(np.hstack((vertices_lores[i,:],1)),np.linalg.inv(xfm_nii).T))
for p,parameter in enumerate(mri_maps):
hippocampus_mri_metrics[i,p] = mri_data[int(idx[0]),int(idx[1]),int(idx[2]),p]
fig = go.Figure(data=[go.Mesh3d(
x=vertices_lores[:,0],
y=vertices_lores[:,1],
z=vertices_lores[:,2],
i=faces_lores[:,0],
j=faces_lores[:,1],
k=faces_lores[:,2],
intensity=hippocampus_mri_metrics[:,0],
colorscale='viridis',
cmax=6,
cmin=15,
showscale=False)])
fig.add_trace(go.Scatter3d(
x=df_nodes.x[components_mask],
y=df_nodes.y[components_mask],
z=df_nodes.z[components_mask],
mode='markers',
marker=dict(
size=diameter_scaled[components_mask],
line_width=0,
opacity=1,
color=node_geometry[components_mask,4],
colorscale='viridis',
cmin=6,
cmax=15,
showscale=True,
colorbar=dict(
title="B1+ (mT)",
titleside="right",
thickness=10
)
)))
fig.update_layout(dict1=vasculature_settings)
fig.write_json("{}/vascular_tree_b1.json".format(vasculature_path))
fig.write_html("{}/vascular_tree_b1.html".format(vasculature_path))
fig.show()
import seaborn as sns
axis_labels = ['Volume (mm$^3$)','Lateral area (mm$^2$)','Surface area (mm$^2$)']
fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10,4))
for g in np.arange(0,3):
sns.scatterplot(x=node_geometry[:,3], y=node_geometry[:,g], alpha=0.25, linewidths=0, ax=axes[g])
sns.despine()
axes[g].set_xlabel('Diameter (mm)')
axes[g].set_ylabel(axis_labels[g])
plt.tight_layout()
# Save network to Panda's pickle
df_nodes.to_pickle('{}/sub-{}_{}_network_nodes.pkl'.format(vasculature_path, subject, hemisphere))
df_edges.to_pickle('{}/sub-{}_{}_network_edges.pkl'.format(vasculature_path, subject, hemisphere))
# Save others to Numpy npy
np.save('{}/sub-{}_{}_network_components.npy'.format(vasculature_path, subject, hemisphere),G_cc)
np.save('{}/sub-{}_{}_hippocampus_vessel.npy'.format(vasculature_path, subject, hemisphere),hippocampus_vessel_metrics_hires)
np.save('{}/sub-{}_{}_hippocampus_mri.npy'.format(vasculature_path, subject, hemisphere),hippocampus_mri_metrics_hires)